home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_5.7 / XSGLL / XSGLL.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  27KB  |  936 lines
  1. /* 
  2.  * xsgll.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * XSGLL
  11.  *
  12.  * construct segm adj lists (SAL) and segm group lists (SGL)
  13.  * by invoking manipulation routines for doubly linked lists;
  14.  * the data item associated with a particular link (node) of a list
  15.  * is defined as struct itemtype (here: Segmtype, see below);
  16.  *
  17.  * one SAL, SAL<i>, is constructed (in the form of a doubly linked 
  18.  * list) for each segment <i> in the image, obtained from a Wall-Damielsson
  19.  * linear approximation and stored in a file of type  .seg, and entered 
  20.  * into an array  sall[] of SALs;
  21.  * for convenience of subsequent construction of SGLs, each SAL 
  22.  * contains segments in the order of decreasing inverse overlap, pji; 
  23.  * the reference element <i> (whose pji is set to 999.0) is placed at 
  24.  * the top; by setting an approriate switch (in sall.c), SALs may also be
  25.  * constructed in the order of increasing dij, (with dij for the
  26.  * reference segment set to 0.0); in that case, each list represents
  27.  * the (forward and backward) segm adjacency list of segment <i>;
  28.  *
  29.  * the segm group list is constructed by contration of SALs
  30.  * and ordered according to dij:
  31.  *   segments j for which dij < 0 will be precede i in the list (b-SAL)
  32.  *   segments j for which dij > 0 will be follow i in the list (f-SAL)
  33.  *
  34.  * note: 
  35.  *   -->short int and int are assumed to be identical (16 bits)!!
  36.  *
  37.  *   -->improvements are possible to substantially reduce the amount
  38.  *   of dynamically allocated memory required for SAL and SGL storage:
  39.  *   in multiple copies of each segment entered in SALs, keep only
  40.  *   the segment index, e.g. <j>, and the pertinent mutual parameters
  41.  *   referring to the reference segment <i>, i.e., dij, pij, pji;
  42.  *   this reduces the amount of memory required per list node from
  43.  *   currently 12+32 bytes to 12+14 bytes;
  44.  *   initialize only non-empty SGLs of which there are only n_sgl<<n_segm;
  45.  *
  46.  * note:
  47.  *   to reduce amount of dynamically allocated memory in the way indicated
  48.  *   in the note above, the present version of xsgll and associated modules 
  49.  *   maintain two segment structures (defined in sgl.h): 
  50.  *   the first ( struct Segm ) contains data relating to only a single segment 
  51.  *   (endpoints, slope, segm. index and line index), while the second structure 
  52.  *   ( struct Segmtype ) serves to store information relating to pairs of 
  53.  *   segments, relevant to SAL and SGL construction
  54.  *   (pij, pji, dij, segm. index, sgl_level, SalStat);
  55.  *
  56.  *   the original and conceptually preferable version of the code maintains only
  57.  *   one structure referring to segments; 
  58.  *   see: lsgl.h, lxsgll.c, 
  59.  *     lsgll.c, lsall.c, ltestsegm.c, lsgl_to_poly.c,lsgl_stats.c
  60.  *
  61.  */
  62.  
  63. #include "xsgll.h"
  64.  
  65. #define MS_TESTSET    13           /* select one of two test sets */
  66. #define LRM_TESTSET    11          /* contained in testsegm.seg */
  67. #define TOTAL_SEG    LRM_TESTSET
  68.  
  69. #undef    DEBUG
  70. #undef    DBG_MEM
  71.  
  72. #undef    ECHO_INPUT
  73. #undef    SHOW_VERT
  74. #undef    SAVE_VERT                /* save cluster vertex coordinates */
  75.  
  76.  
  77.  
  78. /* default values for thresholds used for SAL and SGL construction */
  79. /*
  80.  * segment pairs are processed (in sall.c) if their relative orientation
  81.  * differs by less than ANGL_THRESH, their overlap pij exceeds OVLP_THRESH
  82.  * and their perpendicular distance (see sall.c, prlpar.c) does not
  83.  * exceed D_MAX
  84.  */
  85. #define ANGL_THRESH    10.0        /* threshold in degrees */
  86.        /* tg(th2-th1) is comp. to */
  87.          /* slp_thresh = tg(ANGL_THRESH) */
  88. #define    OVLP_THRESH    50.0        /* normalized to vary between 0.0 amd 100.0 */
  89. #define    D_MAX        50.0             /* in pixels */
  90.  
  91. #define    N_BINS        10              /* number of bins in cluster area histogram */
  92.  
  93.  
  94. #define    NO_DISPLAY    -1
  95. #define    TERSE        0
  96. #define    VERBOSE        1
  97.  
  98. #define    NA_MAX        128             /* max no allowed vertices in ph.c */
  99.  
  100. /*
  101.  * global variables
  102.  */
  103. int loop_switch = 1;            /* enable/disable loop over poly approx */
  104. int hull_switch = 1;            /* enable/disable evaluation of hull eval */
  105. int n_ap_max = NA_MAX;
  106. extern char *optarg;
  107. extern int optind, opterr;
  108.  
  109. int DRAW_MODE = 1;              /* 1: overwrite existing image */
  110.        /* 3: XOR, i. e. overlay output onto img */
  111. int WRITE_FILE = 0;
  112. int DISP_MODE = NO_DISPLAY;     /* VERBOSE, TERSE or NO_DISPLAY */
  113. int TEST_INPUT = 0;
  114. int INPUT = 1;
  115. int COLOR_CODE = 0;             /* when set, segm clusters are color coded */
  116.       /* acc to segm normal dirn (gpoly_hull.c) */
  117. int DISPL_CONV_HULL = 1;        /* when set, display convex hull */
  118. int MERGE = 1;                  /* when set, merge segments after SGL construction */
  119.  
  120.  
  121. /*
  122.  * initialize list parameters
  123.  *
  124.  * note: matching function is set by call to llsetmatch() in code as needed
  125.  */
  126. void
  127. init_sall (struct linklist *segm_adj_list)
  128. {
  129.  
  130.   llzero (segm_adj_list);       /* init empty list */
  131.   llsetsize (sizeof (struct Segmtype), segm_adj_list);
  132. }
  133.  
  134.  
  135. /*
  136.  * initialize list parameters
  137.  *
  138.  * note: matching function is set by call to llsetmatch() in code as needed
  139.  */
  140. void
  141. init_sgll (struct linklist *segm_group_list)
  142. {
  143.  
  144.   llzero (segm_group_list);     /* init empty list */
  145.   llsetsize (sizeof (struct Segmtype), segm_group_list);
  146. }
  147.  
  148. /*
  149.  * deallocate memory occupied by links in a list
  150.  */
  151. int
  152. rm_llistlink (struct linklist *list)
  153. {
  154.   int link_count;
  155.  
  156.   link_count = 0;
  157.   llhead (list);
  158.   do {
  159.     lldelete (list);
  160.     link_count++;
  161.  
  162.   } while (list->listlength != 0);
  163.  
  164.   return (link_count);
  165. }
  166.  
  167.  
  168.  
  169. /*
  170.  * display list (terse version)
  171.  */
  172. void
  173. tshow_segm_list (struct linklist *list, int index)
  174. {
  175.   struct Segmtype *segm;
  176.   int isegm = 0;
  177.  
  178.   printf ("\n...list<%d>\n", index);
  179.   llhead (list);
  180.   if (ll_length (list) != 0) {  /* list not empty */
  181.     do {
  182.       segm = (struct Segmtype *) list->clp->item;
  183.       printf (" %3d", segm->segm_ind);
  184.       if ((isegm + 1) % 12 == 0) {
  185.         isegm = 0;
  186.         printf ("\n");
  187.       }
  188.     } while (llnext (list) == True);
  189.   }
  190.   else
  191.     printf ("  list empty\n");
  192. }
  193.  
  194. /*
  195.  * display list
  196.  */
  197. void
  198. show_segm_list (struct linklist *list, int index)
  199. {
  200.   struct Segmtype *segm;
  201.  
  202.   printf ("\n...list <%d>\n", index);
  203.   llhead (list);
  204.   if (ll_length (list) != 0) {  /* list not empty */
  205.     do {
  206.       segm = (struct Segmtype *) list->clp->item;
  207.       printf ("  segm: %3d", segm->segm_ind);
  208.       printf ("  d(%2d,%2d): %6.2f",
  209.               index, segm->segm_ind, segm->dij);
  210.       printf ("\tp(%2d,%2d): %6.2f",
  211.               index, segm->segm_ind, segm->pij);
  212.       printf ("  p(%2d,%2d): %6.2f",
  213.               segm->segm_ind, index, segm->pji);
  214.       printf ("  SGL lev: %d", segm->sgl_level);
  215.       printf ("\n");
  216.     } while (llnext (list) == True);
  217.   }
  218.   else
  219.     printf ("  list empty\n");
  220. }
  221.  
  222.  
  223.  
  224. /*
  225.  * evaluate slope of linear segment
  226.  */
  227. float
  228. fslope (struct spoint pt1, struct spoint pt2)
  229. {
  230.   double m;
  231.   double max_slope = 500.0;     /* max. poss. slope */
  232.   int x1, y1, x2, y2;
  233.  
  234.   x1 = (int) pt1.x;
  235.   y1 = (int) pt1.y;
  236.   x2 = (int) pt2.x;
  237.   y2 = (int) pt2.y;
  238.  
  239.   if ((int) (x2 - x1) == 0) {
  240.     if ((int) (y2 - y1) == 0)
  241.       m = 0.0;
  242.     else
  243.       m = max_slope;
  244.   }
  245.   else
  246.     m = ((double) (y2 - y1)) / ((double) (x2 - x1));
  247.  
  248.   return ((float) m);
  249. }
  250.  
  251.  
  252.  
  253.  
  254. /*
  255.  * read first line in data file to determine size of record
  256.  */
  257. void
  258. get_record_size (FILE * file, int *n_segm, int *xmax, int *ymax)
  259. {
  260.   int retval;
  261.  
  262.   if (((retval = fscanf (file, "%d %d %d", n_segm, xmax, ymax)) != 3) || (retval == EOF)) {
  263.     printf ("wrong input file format!\n");
  264.     exit (1);
  265.   }
  266. }
  267.  
  268.  
  269. /*
  270.  * read segment data from formatted data file (generally of type .seg)
  271.  * written by pcctoseg.c; first column contains segment index.
  272.  */
  273. void
  274. init_segm (FILE * fp, int n_segm, struct Segm *segm)
  275. {
  276.   int i;
  277.   int retval;
  278.   float skip_slope;
  279.  
  280.  
  281.   for (i = 0; i < n_segm; i++) {
  282.     retval = fscanf (fp, "%4d %3d %3d %3d %3d %f %3d",
  283.                      &((segm + i)->segm_ind),
  284.                      &((segm + i)->ptO.x), &((segm + i)->ptO.y),
  285.                      &((segm + i)->ptF.x), &((segm + i)->ptF.y),
  286.                      &skip_slope, &((segm + i)->line_ind));
  287.  
  288. /* 
  289.  * recompute slope with function internal to this module
  290.  */
  291.     (segm + i)->slope = (float) slope (&(segm + i)->ptO, &(segm + i)->ptF);
  292.  
  293.   }
  294. }
  295.  
  296.  
  297.  
  298.  
  299. /*
  300.  * gprintf() functions:
  301.  * write to std output and, if option WRITE_FILE set, to file wbuf (stream)
  302.  */
  303. void
  304. gprintf (FILE * fpOut, char *fmt,...)
  305. {
  306.   va_list arg_ptr;
  307.  
  308.   va_start (arg_ptr, fmt);
  309.   vprintf (fmt, arg_ptr);
  310.   if (WRITE_FILE == 1)
  311.     vfprintf (fpOut, fmt, arg_ptr);
  312.   va_end (arg_ptr);
  313.  
  314. }
  315.  
  316.  
  317.  
  318. /*
  319.  * error message
  320.  */
  321. void
  322. exitmess (char *prompt, int status)
  323. {
  324.   printf (prompt);
  325.   printf ("\n");
  326.  
  327.   exit (status);
  328. }
  329.  
  330.  
  331. /*
  332.  * usage of routine
  333.  */
  334. void
  335. usage (char *progname)
  336. {
  337.   progname = last_bs (progname);
  338.   printf ("USAGE: %s infile outimg [-t] [-w] [-a ang_thresh][-p ovlp_thresh]\n", progname);
  339.   printf ("                        [-d d_max] [-m] [-o] [-c] [-h] [-L]\n");
  340.   printf ("\n%s constructs segm adj lists (SAL) and segment group lists (SGL)\n", progname);
  341.   printf ("given an input file containing line segment data; writes output\n");
  342.   printf ("data into file and creates image containing line segment clusters\n\n");
  343.   printf ("ARGUMENTS:\n");
  344.   printf ("        infile: input filename (.seg) (ASCII)\n");
  345.   printf ("        outimg: output image filename (TIF)\n\n");
  346.   printf ("OPTIONS:\n");
  347.   printf ("            -t: employ test data in test_segm.c\n");
  348.   printf ("            -w: write sgl parameters and stats to file fn.sgl\n");
  349.   printf ("-a angl_thresh: set angle threshold - default: 10 deg\n");
  350.   printf ("-p ovlp_thresh: set ovlp_thresh (0 <= pij <= 100) - default: 50\n");
  351.   printf ("      -d d_max: set maximum dist between segments - default: 50 pixels\n");
  352.   printf ("            -m: do not merge segments after SGL constr\n");
  353.   printf ("            -o: overlay topological features onto image\n");
  354.   printf ("            -c: display segment clusters\n");
  355.   printf ("            -h: do not display conv hull, only orig. polygon\n");
  356.   printf ("            -L: print Software License for this module\n");
  357.   exit (1);
  358. }
  359.  
  360.  
  361. void
  362. main (int argc, char *argv[])
  363. {
  364.   FILE *file, *fpOut;
  365.   static char *buf =
  366.   {"inpsegm.seg"};              /* default input file name */
  367.   static char *wsuffix =
  368.   {".sgl"};                     /* suffix for output file name */
  369.   static char wbuf[13];
  370.   int ich, is;
  371.  
  372.   /* threshold variables */
  373.   int i_arg;
  374.   float angl_thresh = (float) ANGL_THRESH;
  375.   float ovlp_thresh = (float) OVLP_THRESH;
  376.   float d_max = (float) D_MAX;
  377.  
  378.   /* SAL */
  379.   struct linklist *sall;        /* array of segm adj lists */
  380. #ifdef DEBUG
  381.   struct linklist *csall;       /* ptr to current SAL */
  382. #endif
  383.   struct Segm *inSegm;          /* sample input */
  384.   int n_segm;
  385.   int n_sall = 0;
  386.   int i;
  387.  
  388.   /* SGL */
  389.   struct linklist *sgll;        /* array of segm group lists */
  390.   struct linklist *csgll;       /* ptr to current SGL */
  391.   int n_sgl = 0;
  392.   int rm_link;
  393.  
  394.   /* vertices and convex hull */
  395.   struct spoint *poly_vert, *v_ap;  /* ptr to (cycl) poly */
  396.   int *x_vert, *y_vert;
  397.   int iv, nv, nvpp;
  398.   long lnv;
  399.   struct spoint *hullctr;
  400.   int xhullc, yhullc;
  401.   double av_dirn = (double) 0.0;  /* average directn of SGL segm */
  402.   float s_nem, *clust_dirn;
  403.  
  404.   float maxlen, maxwidth;
  405.   float *asp_ratio;             /* eccentricity of SGL */
  406.   int *ex_hist;
  407.   float min_asp_ratio, max_asp_ratio, ex_bin_width;
  408.  
  409.   int min_cluster_area, max_cluster_area, area_bin_width;
  410.   int ibin, tot_number;
  411.   int *clust_area, *area_hist;
  412.   int *hull_area;
  413.  
  414.   int *clust_lsz, *lsz_hist;
  415.   int min_lsz, max_lsz, nlsz_bins;
  416.  
  417.  
  418.   struct Bdy bd, *bdp = &bd;    /* hsbaird's boundary struct */
  419.   double tot_phi = (double) -8.0;
  420.  
  421.   /* graphics */
  422.   Image *ip;
  423.   int xmax = 0;
  424.   int ymax = 0;
  425.  
  426. /* 
  427.  * cmd line options ( see usage() ):
  428.  */
  429.   static char *optstring = "twa:p:d:mochL";
  430.  
  431.  
  432. /*
  433.  * parse command line
  434.  */
  435.   optind = 3;
  436.   opterr = 1;                   /* give error messages */
  437.  
  438.   if (argc < 3)
  439.     usage (argv[0]);
  440.  
  441.   while ((i_arg = getopt (argc, argv, optstring)) != EOF) {
  442.     switch (i_arg) {
  443.  
  444.     case 't':
  445.       printf ("\n...option %c: ", i_arg);
  446.       printf ("employ test data in test_segm.c\n");
  447.       n_segm = TOTAL_SEG;
  448.       TEST_INPUT = 1;
  449.       INPUT = 0;
  450.       break;
  451.     case 'w':
  452.       printf ("\n...option %c: ", i_arg);
  453.       WRITE_FILE = 1;
  454.       break;
  455.     case 'a':
  456.       printf ("\n...option %c: ", i_arg);
  457.       printf ("set angle threshold to: %f\n",
  458.               angl_thresh = (float) atof (optarg));
  459.       break;
  460.     case 'p':
  461.       printf ("\n...option %c: ", i_arg);
  462.       printf ("set overlap threshold to: %f\n",
  463.               ovlp_thresh = (float) atof (optarg));
  464.       break;
  465.     case 'd':
  466.       printf ("\n...option %c: ", i_arg);
  467.       printf ("set max segm-to-segm dist to: %f\n",
  468.               d_max = (float) atof (optarg));
  469.       break;
  470.     case 'm':
  471.       printf ("\n...option %c:", i_arg);
  472.       printf ("  do not merge segments after SGL construction\n");
  473.       MERGE = 0;
  474.       break;
  475.     case 'o':
  476.       printf ("\n...option %c:", i_arg);
  477.       printf ("  overlay output on img in buf 1\n");
  478.       DRAW_MODE = 3;
  479.       break;
  480.     case 'c':
  481.       printf ("\n...option %c:", i_arg);
  482.       printf ("  display segment clusters\n");
  483.       COLOR_CODE = 1;
  484.       break;
  485.     case 'h':
  486.       printf ("\n...option %c:", i_arg);
  487.       printf ("  do not display convex hull, only original polygon\n");
  488.       DISPL_CONV_HULL = 0;
  489.       break;
  490.     case 'L':
  491.       print_sos_lic ();
  492.       exit (0);
  493.     default:
  494.       printf ("\n...unknown condition encountered\n");
  495.       exit (1);
  496.       break;
  497.     }
  498.   }
  499.  
  500.  
  501.   if (!TEST_INPUT) {
  502.     printf ("read segment data from file %s\n", buf = argv[1]);
  503.     /* get size of data set from file */
  504.     if ((file = fopen (buf, "r")) == NULL) {
  505.       printf ("\n...could not open file %s\n", buf);
  506.       exit (1);
  507.     }
  508.     /* initialize size variables */
  509.     get_record_size (file, &n_segm, &xmax, &ymax);
  510.     if (n_segm > MAX_REC_SIZE) {
  511.       printf ("\n...record size of %d exceeds limit of %d\n", n_segm, MAX_REC_SIZE);
  512.       exit (1);
  513.     }
  514.  
  515.     /* construct output file name */
  516.     ich = 0;
  517.     while ((*(wbuf + ich) = *(buf + ich)) != '.')
  518.       ich++;
  519.     for (is = 0; is < 4; is++)
  520.       *(wbuf + (ich) + is) = *(wsuffix + is);
  521.  
  522.     if (WRITE_FILE) {
  523.       if ((fpOut = fopen (wbuf, "w")) == NULL) {
  524.         printf ("\n...cannot open %s for writing\n", wbuf);
  525.         exit (1);
  526.       }
  527.       gprintf (fpOut, "logging SGL params in file %s\n", wbuf);
  528.     }
  529.   }
  530.  
  531.  
  532. #ifdef DBG_MEM
  533.   gmod_logfile (1);             /* turn on(1)/off(0) gmod log file */
  534. #endif
  535.  
  536. /*
  537.  * display parameter settings employing global printf function;
  538.  * open log file
  539.  */
  540.   gprintf (fpOut, "\nthreshold settings:\n");
  541.   gprintf (fpOut, "    ANGL_THRESH = %f  slp_thresh = %f\n",
  542.            angl_thresh, tan (angl_thresh * (2.0 * PI / 360.0)));
  543.   gprintf (fpOut, "    OVLP_THRESH = %f\n", ovlp_thresh);
  544.   gprintf (fpOut, "    DIST_THRESH = %f\n", d_max);
  545.  
  546.  
  547. /*
  548.  * memory allocation
  549.  */
  550.   if ((inSegm = (struct Segm *) calloc (n_segm, sizeof (struct Segm))) == NULL)
  551.       exitmess ("can't allocate mem for struct Segm *inSegm", 1);
  552.   if ((sall = (struct linklist *) calloc (n_segm, sizeof (struct linklist))) == NULL)
  553.       exitmess ("can't allocate mem for struct linklist *sall", 1);
  554.   if ((sgll = (struct linklist *) calloc (n_segm, sizeof (struct linklist))) == NULL)
  555.       exitmess ("can't allocate mem for struct linklist *sgll", 1);
  556.  
  557.  
  558. /*
  559.  * fetch segment data
  560.  */
  561.   if (TEST_INPUT == 1) {
  562.     printf ("\n...initialize test segments\n");
  563.     init_testsegm (inSegm, n_segm);
  564.     ymax = HEIGHT;
  565.     xmax = WIDTH;
  566.   }
  567.   else if (INPUT == 1) {
  568.     gprintf (fpOut, "input file for segment data: %s\n", buf);
  569.     init_segm (file, n_segm, inSegm);
  570.     fclose (file);
  571.   }
  572.  
  573. /*
  574.  * initialize graphics
  575.  */
  576.   ip = ImageAlloc ((long) ymax, (long) xmax, BPS);
  577.  
  578.  
  579. #ifdef ECHO_INPUT
  580.   printf ("\n...segment data, asread in from file:\n");
  581.   for (i = 0; i < n_segm; i++) {
  582.     printf ("%4d  %3d %3d %3d %3d\n",
  583.             (inSegm + i)->segm_ind,
  584.             (inSegm + i)->ptO.x, (inSegm + i)->ptO.y,
  585.             (inSegm + i)->ptF.x, (inSegm + i)->ptF.y);
  586.   }
  587. #endif
  588.  
  589.   printf ("\n...initialize %d empty Segm Adj Lists (SAL)\n", n_segm);
  590.   for (i = 0; i < n_segm; i++)
  591.     init_sall (sall + i);
  592.  
  593. #ifdef DEBUG
  594.   csall = sall + 0;
  595.   printf ("\n...structure parameters for segm adj list number 0:\n");
  596.   printf ("sizeof(struct linklist): %d\n", sizeof (struct linklist));
  597.   printf ("        item length: %d\n", csall->itemlength);
  598.   printf ("    cur list length: %d\n", csall->listlength);
  599. #endif
  600.  
  601.  
  602.   printf ("...construct SAL for each of %d segments\n", n_segm);
  603.   printf ("     (in the form of a doubly linked list of segments)\n");
  604.   n_sall = construct_sall (sall, inSegm, n_segm,
  605.                            angl_thresh, ovlp_thresh, d_max);
  606.  
  607. /*
  608.  * display results
  609.  */
  610.   printf ("\n...done with construction of pji-ordered SAL:\n");
  611.  
  612.  
  613.   if (DISP_MODE == VERBOSE)
  614.     for (i = 0; i < n_segm; i++)
  615.       show_segm_list (sall + i, i);
  616.   else if (DISP_MODE == TERSE)
  617.     for (i = 0; i < n_segm; i++)
  618.       tshow_segm_list (sall + i, i);
  619.  
  620.  
  621. /*
  622.  * SGL construction
  623.  */
  624.   printf ("\n...initialize %d empty Segm Group Lists (SGL)\n", n_segm);
  625.   for (i = 0; i < n_segm; i++)
  626.     init_sgll (sgll + i);
  627.  
  628.  
  629. #ifdef DEBUG
  630.   csgll = sgll + 0;
  631.   printf ("\n...structure parameters for segm adj list number 0:\n");
  632.   printf ("        item length: %d\n", csgll->itemlength);
  633.   printf ("    cur list length: %d\n", csgll->listlength);
  634. #endif
  635.  
  636.   n_sgl = construct_sgll (inSegm, sall, sgll, n_segm, (double) ovlp_thresh);
  637.  
  638. /*
  639.  * display results
  640.  */
  641.   printf ("\n...done with construction of dij-ordered SGL:\n");
  642.  
  643.   gprintf (fpOut, "\tnumber of SGLs: %3d\n", n_sgl);
  644.   gprintf (fpOut, "\t          ACCEPT_LEVEL: %3d\n", ACCEPT_LEVEL);
  645.   gprintf (fpOut, "\t          min sz for proc: %3d\n", MIN_SGL_SIZE);
  646.  
  647.   if (DISP_MODE == VERBOSE)
  648.     for (i = 0; i < n_segm; i++)
  649.       show_segm_list (sgll + i, i);
  650.   else if (DISP_MODE == TERSE)
  651.     for (i = 0; i < n_segm; i++)
  652.       tshow_segm_list (sgll + i, i);
  653.  
  654.  
  655. /*
  656.  * construct an ordered, cyclic sequence of vertices from the
  657.  * endpoints of the segments in a segm group list; handle each
  658.  * SGL in turn, skipping those containing less than MIN_SGL_SIZE segments;
  659.  *
  660.  * there will be n_sgl <= n_sall <= n_segm (closed) polygons, each containing 
  661.  * nv+1 points: nv thus represents the number of distinct vertices!
  662.  */
  663. #ifdef DEBUG
  664.   printf ("\n...size of struct Bdy: %d\n", sizeof (struct Bdy));
  665. #endif
  666.  
  667.   if ((clust_dirn = (float *) calloc (n_segm, sizeof (float))) == NULL)
  668.       exitmess ("can't allocate mem for clust_dirn", 1);
  669.  
  670.   if ((clust_area = (int *) calloc (n_segm, sizeof (int))) == NULL)
  671.       exitmess ("can't allocate mem for clust_area", 1);
  672.  
  673.   if ((clust_lsz = (int *) calloc (n_segm, sizeof (int))) == NULL)
  674.       exitmess ("can't allocate mem for clust_lsz", 1);
  675.  
  676.   if ((hull_area = (int *) calloc (n_segm, sizeof (int))) == NULL)
  677.       exitmess ("can't allocate mem for hull_area", 1);
  678.  
  679.   if ((asp_ratio = (float *) calloc (n_segm, sizeof (float))) == NULL)
  680.       exitmess ("can't allocate mem for asp_ratio", 1);
  681.  
  682.  
  683.  
  684.   min_cluster_area = 10000;
  685.   max_cluster_area = 0;
  686.   min_asp_ratio = (float) 1000.0;
  687.   max_asp_ratio = (float) 0.0;
  688.   min_lsz = 1000;
  689.   max_lsz = 0;
  690.   n_sgl = 0;
  691.   for (i = 0; i < n_segm; i++)
  692.     *(clust_dirn + i) = (float) 999.0;
  693.  
  694.   for (i = 0; i < n_segm; i++) {  /* loop over SGLs */
  695.     csgll = sgll + i;
  696.     if ((nv = 2 * csgll->listlength) >= 2 * MIN_SGL_SIZE) {
  697.       n_sgl++;                  /* count SGLs of sufficient length */
  698.  
  699.       if ((poly_vert = (struct spoint *) calloc (nv + 1, sizeof (struct spoint))) == NULL)
  700.           exitmess ("can't allocate mem for poly_vert", 1);
  701.  
  702.       if ((v_ap = (struct spoint *) calloc (nv + 1, sizeof (struct spoint))) == NULL)
  703.           exitmess ("\n...mem allocation for v_ap failed\n", 1);
  704.  
  705.       if ((x_vert = (int *) calloc (nv + 1, sizeof (int))) == NULL)
  706.           exitmess ("can't allocate mem for x_vert", 1);
  707.       if ((y_vert = (int *) calloc (nv + 1, sizeof (int))) == NULL)
  708.           exitmess ("can't allocate mem for y_vert", 1);
  709.  
  710.  
  711. /* 
  712.  * construct polygon (closed array of nv+1 vertices) and 
  713.  * determine average segment orientation; note that merging
  714.  * of segments after SGL construction will result in the elimination
  715.  * of segments; for each segment deleted from the SGL, the number of 
  716.  * vertices must be diminished by two
  717.  */
  718.       nvpp = nv + 1;
  719.       av_dirn = construct_va (inSegm, csgll, &nvpp, poly_vert);
  720.       *(clust_dirn + i) = (float) av_dirn;
  721.       nv = nvpp - 1;
  722.  
  723. /* 
  724.  * evaluate cluster longitudinal size (no of segments) and aspect ratio, 
  725.  * a cluster being defined as the set of linear segments in a SGL
  726.  */
  727.       *(clust_lsz + i) = nv / 2;
  728.       if (*(clust_lsz + i) > max_lsz)
  729.         max_lsz = *(clust_lsz + i);
  730.       if (*(clust_lsz + i) < min_lsz)
  731.         min_lsz = *(clust_lsz + i);
  732.  
  733.       *(asp_ratio + i) = eccentr (inSegm, csgll,
  734.                                   &maxlen, &maxwidth);
  735.       if (*(asp_ratio + i) > max_asp_ratio)
  736.         max_asp_ratio = *(asp_ratio + i);
  737.       if (*(asp_ratio + i) < min_asp_ratio)
  738.         min_asp_ratio = *(asp_ratio + i);
  739.  
  740.  
  741. /* 
  742.  * collect cluster areas to be used in histogram evaluation and 
  743.  * determine minimum and maximum cluster areas
  744.  */
  745.       for (iv = 0; iv <= nv; iv++) {
  746.         *(x_vert + iv) = (poly_vert + iv)->x;
  747.         *(y_vert + iv) = (poly_vert + iv)->y;
  748.       }
  749.  
  750.       *(clust_area + i) = zero_moment (nv, x_vert, y_vert);
  751.       if (*(clust_area + i) > max_cluster_area)
  752.         max_cluster_area = *(clust_area + i);
  753.       if (*(clust_area + i) < min_cluster_area)
  754.         min_cluster_area = *(clust_area + i);
  755.  
  756. /*
  757.  * display geometrical parameters
  758.  */
  759. #ifdef DEBUG
  760.       gprintf (fpOut, "\ngeometr prop of cluster %d:\n", i);
  761.       gprintf (fpOut, "   no of segments: %d\n",
  762.                csgll->listlength);
  763.       gprintf (fpOut, "   area = %d\n", *(clust_area + i));
  764.       gprintf (fpOut, "   length = %6.2f, width = %6.2f\n",
  765.                maxlen, maxwidth);
  766.       gprintf (fpOut, "   asp_ratio = %6.2f\n", *(asp_ratio + i));
  767.       gprintf (fpOut, "   avrg_dirn = %6.2f\n", av_dirn);
  768. #endif
  769.  
  770.  
  771. #ifdef SHOW_VERT
  772.       printf ("\n...ordered vertex array for SGL<%d>\n", i);
  773.       printf ("    n_vertex: %d\n", nv);
  774.       printf ("      avrg slope: %lf\n", av_dirn);
  775.       for (iv = 0; iv <= nv; iv++)
  776.         printf ("    (%d, %d)\n",
  777.                 (poly_vert + iv)->x, (poly_vert + iv)->y);
  778.  
  779. #endif
  780. #ifdef SAVE_VERT
  781.       gprintf (fpOut, "\n...ordered vertex array for SGL<%d>\n", i);
  782.       gprintf (fpOut, "    n_vertex: %d\n", nv);
  783.       gprintf (fpOut, "      avrg slope: %lf\n", av_dirn);
  784.       for (iv = 0; iv <= nv; iv++)
  785.         gprintf (fpOut, "    (%d, %d)\n",
  786.                  (poly_vert + iv)->x, (poly_vert + iv)->y);
  787. #endif
  788.  
  789.  
  790.  
  791. /* 
  792.  * construct convex hull of cluster
  793.  */
  794.       lnv = (long) nv;
  795.       fill_bdy_structure (bdp, poly_vert, lnv, tot_phi);
  796.       hullctr = poly_hull_tol (bdp, v_ap, av_dirn, hull_area + i, (float) 0.0, ip, WHITE);
  797.  
  798.       xhullc = hullctr->x;
  799.       yhullc = hullctr->y;
  800.  
  801. /*
  802.  * deallocate memory
  803.  */
  804.       free (x_vert);
  805.       free (y_vert);
  806.       free (poly_vert);
  807.       free (v_ap);
  808.  
  809.       free (bdp->hpp);          /* allocated in poly_approx.c */
  810.       free (bdp->ap);           /* allocated in poly_approx.c */
  811.       free (bdp->v);            /* allocated in poly_edge.c */
  812.     }
  813.   }
  814.  
  815.   if (!n_sgl) {
  816.     printf ("No clusters found to process\nExiting!\n");
  817.     exit (0);
  818.   }
  819.  
  820.  
  821.  
  822. /*
  823.  * evaluate SGL (cluster) longitudinal size (no of segments) histogram
  824.  */
  825.   nlsz_bins = max_lsz - min_lsz + 1;
  826.   if ((lsz_hist = (int *) calloc (nlsz_bins, sizeof (int))) == NULL)
  827.       exitmess ("can't allocate mem for lsz_hist", 1);
  828.   find_lsz_hist (n_segm, clust_lsz, lsz_hist, nlsz_bins, min_lsz);
  829.  
  830.   tot_number = 0;
  831.   gprintf (fpOut, "\nhistogram of cluster longit. sizes:\n");
  832.   gprintf (fpOut, "  maximum no of segm: %d\n", max_lsz);
  833.   gprintf (fpOut, "  minimum no of segm: %d\n", min_lsz);
  834.   for (ibin = 0; ibin < nlsz_bins; ibin++) {
  835.     gprintf (fpOut, " %3d ", *(lsz_hist + ibin));
  836.     if ((ibin + 1) % 10 == 0)
  837.       gprintf (fpOut, "\n");
  838.  
  839.     tot_number += *(lsz_hist + ibin);
  840.   }
  841.   gprintf (fpOut, "\n");
  842.  
  843. /*
  844.  * evaluate SGL (cluster) aspect ratio (eccentricity: ex ) histogram
  845.  */
  846.   ex_bin_width = (max_asp_ratio - min_asp_ratio) / (float) N_BINS;
  847.   if ((ex_hist = (int *) calloc (N_BINS, sizeof (int))) == NULL)
  848.       exitmess ("can't allocate mem for ex_hist", 1);
  849.   find_ex_hist (n_segm, asp_ratio, ex_hist, ex_bin_width, min_asp_ratio);
  850.  
  851.   tot_number = 0;
  852.   gprintf (fpOut, "\nhistogram of cluster aspect ratios:\n");
  853.   gprintf (fpOut, "   maximum asp ratio: %f\n", max_asp_ratio);
  854.   gprintf (fpOut, "   minimum asp ratio: %f\n", min_asp_ratio);
  855.   gprintf (fpOut, "        ex_bin_width: %f\n", ex_bin_width);
  856.   for (ibin = 0; ibin < N_BINS; ibin++) {
  857.     gprintf (fpOut, " %3d ", *(ex_hist + ibin));
  858.     tot_number += *(ex_hist + ibin);
  859.   }
  860.   gprintf (fpOut, "\n");
  861.  
  862. /*
  863.  * evaluate SGL (cluster) area histogram
  864.  */
  865.   area_bin_width = (max_cluster_area - min_cluster_area) / N_BINS;
  866.   if ((area_hist = (int *) calloc (N_BINS, sizeof (int))) == NULL)
  867.       exitmess ("can't allocate mem for area_hist", 1);
  868.   find_area_hist (n_segm, clust_area, area_hist,
  869.                   area_bin_width, min_cluster_area);
  870.  
  871.   tot_number = 0;
  872.   gprintf (fpOut, "\nhistogram of cluster areas:\n");
  873.   gprintf (fpOut, "    maximum area: %d\n", max_cluster_area);
  874.   gprintf (fpOut, "    minimum area: %d\n", min_cluster_area);
  875.   gprintf (fpOut, "      area_bin_width: %d\n", area_bin_width);
  876.   for (ibin = 0; ibin < N_BINS; ibin++) {
  877.     gprintf (fpOut, " %3d ", *(area_hist + ibin));
  878.     tot_number += *(area_hist + ibin);
  879.   }
  880.   gprintf (fpOut, "\n");
  881.  
  882. /*
  883.  * evaluate nematic order parameter for ensemble of clusters
  884.  */
  885.   s_nem = nematic_op (n_sgl, n_segm, clust_dirn);
  886.   gprintf (fpOut, "\nnematic order parameter: S = %6.2f\n", s_nem);
  887.  
  888.  
  889.   free (asp_ratio);
  890.   free (ex_hist);
  891.   free (clust_area);
  892.   free (area_hist);
  893.   free (clust_lsz);
  894.   free (lsz_hist);
  895.   free (clust_dirn);
  896.  
  897.   free (hull_area);
  898. /*
  899.  * deallocate memory assigned to lists
  900.  */
  901.   printf ("\n...deallocating memory for linked lists\n");
  902.   printf ("    SALs:\n");
  903.   for (i = 0; i < n_segm; i++) {
  904.     if ((sall + i)->listlength != 0) {
  905.       rm_link = rm_llistlink (sall + i);
  906. #ifdef DBG_MEM
  907.       printf ("   for SAL<%d>: %d links removed\n",
  908.               i, rm_link);
  909. #endif
  910.     }
  911.   }
  912.   printf ("    SGLs\n");
  913.   for (i = 0; i < n_segm; i++) {
  914.     if ((sgll + i)->listlength != 0) {
  915.       rm_link = rm_llistlink (sgll + i);
  916. #ifdef DBG_MEM
  917.       printf ("   for SGL<%d>: %d links removed\n",
  918.               i, rm_link);
  919. #endif
  920.     }
  921.   }
  922.  
  923.   free (inSegm);                /* free array of input segments */
  924.   free (sall);
  925.   free (sgll);
  926.  
  927.   if (WRITE_FILE == 1)
  928.     fclose (fpOut);
  929.   ImageOut (argv[2], ip);
  930.  
  931. #ifdef DBG_MEM
  932.   printf ("\n...no of bytes still in use: %d\n", gmod_inuse ());
  933.   gmod_logfile (0);
  934. #endif
  935. }
  936.